PI-CONTROLLED BIOREACTOR AS A 
GENERALIZED LIENARD SYSTEM 



V. Ibarra- Junquera and H.C. Rosu 

Potosinian Institute of Science and Technology (IPICyT), 
Apdo Postal 3-74 Tangamanga, 78231 San Luis Potosi, Mexico 



Abstract 

It is shown that periodic orbits can emerge in Cholette's bioreactor model working 
under the influence of a Pi-controller. We find a diffeomorphic coordinate trans- 
formation that turns this controlled enzymatic reaction system into a general- 
ized Lienard form. Furthermore, we give sufficient conditions for the existence and 
uniqueness of limit cycles in the new coordinates. We also perform numerical simu- 
lations illustrating the possibility of the existence of a local center (period annulus). 
A result with possible practical applications is that the oscillation frequency is a 
function of the integral control gain parameter. 
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1 Introduction 



Mixing, understood as interpenetration of particles in different zones of a 
given volume, is an important natural as well as technological process. This 
is even more so when biochemical reactions get involved. For the case of con- 
tinuous stirred tank reactors (CSTRs), Lo and Cholette [1] developed a non- 
ideal isothermal mixing model using a Haldane type chemical reaction rate 
(which is similar to the Monod function for low concentrations but includes 
the inhibitory effect at high concentrations). This model has been studied ex- 
tensively later by many authors ([2], [3], [4], [5], [6]). In particular, Sree and 
Chidambaram [4], [5] focused on the control problem by means of a propor- 
tional integral (PI) control for this case. Indeed, the PI controller is broadly 
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used in the chemical and biochemical industry. Therefore its closed-loop be- 
havior is of much interest. In this paper, we present a novel mathematical 
feature of this closed-loop enzymatic reaction system, namely the possibility 
to be represented as a dynamical system corresponding to a non polynomial 
Lienard oscillator. That means that given a Pi-controlled CSTR governed by 
the usual two-dimensional smooth dynamical system 



we are able to find a diffeomorphic coordinate transformation (Eq. 7 below) 
that allows us to put it into the well-known generalized Lienard form 



where g(X) is continuous on an open interval (a±, bi), the functions F(X) and 
4>(Y) are continuously differentiable on the open intervals (a l5 bi) and (02, b 2 ), 
respectively. In fact, these intervals can be extended to —00 < <2j < < bi < 
00, i = 1, 2. 

In this paper, we show that the Pi-controlled Cholette's CSTR model belongs 
to this class of generalized Lienard systems. Once doing this, we make use of 
the beautiful results encountered in this research area to study the periodic 
solutions near stationary points for this particular application. Besides, the 
Hopf bifurcation is an efficient way to study the existence of periodic orbits. 
In this case, a pair of complex eigenvalues is assumed to exist and to cross 
transversally the imaginary axis. Nevertheless, the fact that a Hopf bifurcation 
guarantees the existence of a limit cycle does not imply its uniqueness [7]. It 
is here where the uniqueness result for Lienard systems comes into play. 

Thus, we extend the class of Lienard-type system to the interesting case of 
Pi-controlled bioreactors for which the results on the existence of limit cycles 
and their number as a function of the control variables could be exploited in 
industrial applications. In general, the study of oscillatory behavior in biore- 
actors is a very important issue since it is generated by the coupled dynamics 
of the most popular controller in industry (the PI one) and the kinetics of 
the biochemical reactions. In addition, we shed light here on an explicit ex- 
ample of a closed-loop system which is of Lienard-type. Since the most direct 
way to interact with the Pi-controlled Cholette's CSTR is through the gain of 
the controller the present analysis provide the users with definite conditions 
for inducing oscillatory behaviors, which is instructive from the pedagogical 
standpoint as well. 

The paper is organized as follows. In Section 2, we discuss the Pi-controlled 



X = P(X,Y), Y = Q(X,Y) 



(1) 



X = <f>(Y)-F(X) 
Y = -g(X), 



(2) 
(3) 
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Cholette's bioreactor model and its basic assumptions. In Section 3, we present 
the coordinate transformation that leads to the Lienard representation of this 
type bioreactor. The existence of limit cycles is discussed in Section 4, and the 
uniqueness consideration are included in Section 5. The numerical simulation 
that we performed indicating the presence of the period annulus are shortly 
described in Section 6. Finally, we end up the work with several concluding 
remarks. 



2 Cholette's Dynamical Model 



The dynamical behavior without control actions (i.e., open-loop operation) 
is governed by a unique nonlinear ordinary differential equation (see Eq. 4). 
The non ideal mixing can by described by the Cholette model [2] . This model 
was studied by Chidambaram [4], who proposed a tuning method for a PI- 
controller. Examples, where this kind of kinetics occurs, can be found in [6]. 
The reactor model is given by the following equations 



nF 



mV 



,, (4) 



where the meaning of the parameters are given in table 1. 

The following assumptions hold in Eq. (4): (i) all model parameters and physic- 
ochemical properties are constant (ii) the reaction occurs in an nonideal mixed 
CSTR, operated under isothermal conditions. The fraction of the reactant feed 
that enters the region of perfect mixing is denoted by n, whereas m denotes 
the fraction of the total volume of the reactor where perfect mixing is achieved. 
For m and n both equal to 1, the system is ideally mixed. The values of the 
parameters m and n can be obtained from the residence time distribution [2]. 
Fig. 1 shows the schematic diagram of the bioreactor configuration modeled 
by Eq. (4). 

Following the previous authors, we consider Q as the manipulated variable (i.e. 
C/ = u) and let ( be the controlled variable [4] . We are especially interested in 
the induced oscillatory behavior of the bioreactor. The common control law in 
this case is of the proportional integral type, that requires a dynamical error 
extension in order to build the closed-loop two dimensional system. Thus, the 
control law is given by 



u± I —K. 



Error — Ki J Error dt^j 
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Table 1 

Variables and parameters of Cholette's model. 



Symbol 


Meaning 


Units 


XT' 

X 


Substrate concentration 


[K mol/m°\ 


I 


Integrated error 


[K mol s/m \ 


771 

b 


Feed flow rate 


[m 3 /s] 


T 7 

V 


Volume 


[m 3 ] 




OF 


Substrate feed concentration 


[K mol fm \ 




Maximal kinetic rate 


ri /„i 
ll/*J 


J - x Z 


Inhibition parameter 


[m? j ' Kniol] 


n 


Mixing parameter 


[dimensionless] 


rn 


Mixing parameter 


[dimensionless] 


Kc 


Proportional gain controller 


[dimensionless] 


K { 


Integral gain controller 


[dimensionless] 


u 


Control input 


[dimensionless] 



Cf Cf 

F 



111 1 



mV 



rruvw 



A 

k 



Cf 

Fig. 1. Schematics of a classical continuous stirred tank bioreactor with imperfect 
mixing corresponding to Cholette's model. 

where K c and Ki are the control gain values. For the sake of simplicity, the 
closed-loop system is written as: 



X 



-X c + 



Ki 

1 + K 2 xy 



C {-K C {X - Ref) - KiY) 



Y = X-Ref 



(5) 
(6) 
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where 



C 



n F 
m V 



x = C 



Eq. (5) describes the dynamical behavior of the concentration, while Eq. (6) 
refers to the dynamical behavior of the integrated error. 



3 Transformation to the Lienard form 



In this section, we show that the system given by Eqs. (5-6) can be rewritten 
as a system of the form (2)- (3), i.e., in the Lienard generalized form. This is 
one of the main results of this work. 

Proposition 1. Under the transformation 



*(X,Y) 



X — x p 

-Y + Y p 



(7) 



where 



Y P = - 



Ref (C + K x + Ref CK 2 (K 2 Ref + 2)) 



CKi (1 + K 2 RefY 



X p = Ref, 



system given by Eq. (5-6) can be written in the generalized Lienard form. With the 
following additional properties 

[Al] 0(0) = and xg(x) > for x / 0; 

[A2] 0(0) = and <jf(x) > for a 2 < y < b 2 ; 

[A3] The curve 4>(y) = F(x) is well defined for all x G (a±, b\). 

Proof. If we substitute X = x + X p and Y = — y + Y p in Eqs. (5) and (6), and we 
choose 



F(x)=x \ C (1 + K c ) + K x 

<f*(y) = CK iy 
g(x)=x , 



1 _ K 2 Ref(2+K 2 (x+2Ref)) \ 



{l+K 2 RefY 



(l + K 2 (x + Ref)Y 



(8) 

(9) 
(10) 
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we get the generalized Lienard form of the Pi-controlled Cholette system. The 
properties [Al], [A2] and [A3] are straightforwardly checked in Eqs. (8)-(10). 

For ty(X, Y) to be a diffeomorphism on the region f2, it is necessary and 
sufficient that the Jacobian d^(X, Y) be nonsingular on Q, and moreover 
that ty(X,Y) is one to one from Q to Since ty(X,Y) is linear, it is 

one to one and the determinant of the Jacobian matrix is constant, then is 
nonsingular in the region Q = [—00,00] x [—00, 00]. 

In the literature, the properties [Ai] are standard properties assumed for 
Lienard systems [10]. We point out that the huge existing literature on Lienard 
systems deals mainly with cases in which F(x) is polynomial [14], [12], [13], 
whereas we are in a case in which F(x) is a nonlinear rational function. Such 
cases are far less studied and there are still many open problems. 



4 Existence of Limit Cycles 



We briefly recall some basic results of the theory of bifurcations of vector 
fields. Roughly speaking, a bifurcation is a change in equilibrium points, peri- 
odic orbits, or in their stability properties, when varying a parameter known 
as bifurcation parameter. The values of the parameter at which the changes 
occur are called bifurcation points. A Hopf bifurcation is characterized by a 
pair of complex conjugate eigenvalues crossing the imaginary axis. Now, sup- 
pose that the dynamical system X = f(X,fj) with X e W 1 and /x e K. has 
an equilibrium point at X eq , for some /i = fi H ; that is / (x eq , /j H ^j = 0. 

Let A(n) = — g X ' - be the Jacobian matrix of the system at the equi- 
librium point. Assume that A has as single pair of purely imaginary 

eigenvalues S (a^) = with oj h > and that these eigenvalues are the 

only ones with the properties ^St(S) = 0. If the following condition is fulfilled 
/ q tne jj pf bifurcation theorem states that at least one limit 

cycle is generated at (x eq , (see [8]). The condition (4) is known as the 
transversality hypothesis. Considering now the nth degree characteristic poly- 
nomial X(S) = S n + aiS*™ -1 + a2<5 , ™~ 2 + . . - + a n , where all the real Oj coefficients 
are positive allows the construction of a Hurwitz matrix Ti, nxn Then one has 
the basic result that the characteristic polynomial is stable if and only if the 
leading principal minors of 7i n xn are all positive [9]. 

To search the Hopf bifurcation we have to calculate the equilibrium points of 
the system (2)-(3) making equal to zero the right-hand-side of the equation, 
and taking as bifurcation parameters the values of the Pi-control gains. Then, 
finding the solution with respect to the state vector x we notice that the 
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closed-loop system has a unique equilibrium point located at the origin. 
Proposition 2. If the parameter K c is such that 



K?-2^<K e <K? , 

then the Lienard system (2) and (3) with the functions given in Eqs. (8)-(10) has 
at least one limit cycle. The upper limit is defined in Eq. (15). 

Proof. If we use the Hurwitz criterion to guarantee that the unique equilibrium 
point is unstable, we evaluate the Jacobian of the system at the origin 



J(0,K c ) 



C ~ (l+K 2 Reff + 2 (l+K 2 Reff ~ CKc C ^ 



-1 







(11) 



Then, from the determinant IS — J(0,K c ), we get the characteristic polynomial 
X(S) = S 2 + aiS + a>2, where 



(l + K 2 Reff V ; 

a 2 = CKi . (13) 

Then, the Hurwitz matrix is given by 



n = 



ai o 
1 a 2 



(14) 



and its principal minors are Hi (K c ) = a\ and H2 (K c ) = a±a2- From this formulas, 
we can note that the stability of the unique equilibrium point depends on the sign 
of Eq. (12). Since all the parameters in the Eqs. (12) and (13) are positive, we can 
induce the local stability as a function of the values of the controller gains involved 
in these equations and then we obtain the bifurcation point as the trivial solution 
of Eq. (12) for K c and restricting Ki > 



C{l + K 2 Reff 

In order to test the transversality condition, the behavior of the eigenvalues of 
i7(0, K c ) in the neighborhood of should be analyzed. Thus, we take K c as 
+ e, with e G R. The transversality condition will be fulfilled if the sign of the 
equations TLi and TL2 changes when the sign of e changes. Substitution of Eq. (15) in 
the principal minor expressions gives Hi (K^) = e (C) and H2 {K^) = e (KiC 2 ) . 
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From the above equations we can appreciate that if we want to have a positive real 
part of the eigenvalues of the Jacobian matrix (11), we need that e be negative. 
In other words, it is below the value given by Eq. (15) where the limit cycles are 
generated. Using , the eigenvalues of the matrix (11) are given by the roots of 
the characteristic polynomial A(5), which are 



eC Ve*C*-4CKj 

51 = ~^ + 2 (16) 

eC Ve*C*-4CKj 

52 = ~~2 2 ' (17) 

where we can notice that the eigenvalues are complex with positive real parts, when 
> e > -2 and K t > 0. 

Eq. (15) is of main importance, because it corresponds to the Hopf bifurcation 
and therefore lies in a neighborhood of the value of the parameter where at 
least one limit cycle is generated. Note, that the Hopf bifurcation by itself can 
not guarantee the uniqueness of the limit cycle, because more than one limits 
cycle could appear [7]. Then we need to use additional constraints in order to 
find the condition for uniqueness. 



5 Uniqueness of Limit Cycles 



Xiao and Zhang [10] gave an interesting theorem on the uniqueness of limit 
cycles for generalized Lienard systems, under the conditions [Al], [A2] and [A3] 
that allows us to prove a novel property of Pi-controlled Cholette bioreactors. 

Theorem 1. Using the notations G(x) = J'q g{x)dx and fix) = F'(x), suppose 
that the system (2) - (3) satisfies the following conditions: 

(i) there exist x\ and x 2 , a\ < x 2 < < x\ < b\ such that F(x\) = F(0) = 0, 
F(x 2 ) > and G(x 2 ) < G(xi); xF(x) < for x 2 < x < x 1: F'(x) > for 
a± < x < x 2 or xi < x < b±, and F(x) 7^ for <| x |<S 1. 

(ii) F(x)f(x)/g(x) is nondecreasing for x\ < x < b±. 
(iii) (j)'(y) is nonincreasing as | y \ increases. 

Then the system (2) and (3) has at most one limit cycle, and it is stable if it exists. 

Note that the above result does not guarantee the existence of limit cycles by 
itself. 

Proposition 3. System (2) and (3) with the functions given by the Eqs. (8)-(10) 
and with Ref > 2/K 2 and K c = K C H + e*, where 
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e = _K l{ -2 + K 2 Ref) 
2C(l + K 2 Ref) 3 

has a unique limit cycle. 

Proof. The existence is given by the Proposition 2, and the uniqueness will be 
proved by showing that the system fulfills the conditions of Theorem 1. Let us take 
Ref* = 2k/K 2 with k > 1 which agrees with the condition stated in Proposition 3. 
Then 



K c H (e*,Ren = -l+ Kl * 



C(l + 2k) 3 ' 

For these values of the parameters, the function F(x) is simplified to the form 



xKx (3 k + k K 2 2 x 2 -4 k 3 + 1) 
(1 + 2k) 3 (1 + K 2 x + 2kY 



The real roots of F*(x) are 0, x\, and —x±, where 
xi = y^FE5(l±2K) 

K^2 



which fulfill the property F{x\) = F(0) = 0. Moreover, taking x 2 = —£xi, where 

< £ < 1, one can see that x 2 < < x\. To evaluate F*(x) in the intervals 
x 2 < x < and < x < x\, it is sufficient to substitute in Eq. (19) x = (3xi, where 

1 = 1, 2 and < /? < 1. Then we can easily verify that 



0>F . ( ^ l)= (zl±^^lkWEB+3 (21) 



(1 + 2k) 2 (k + /V^ (-1 + k)) K 2 



(-1 + k) (0 2 -1) KK x aJ~k (-1 + k) , . 

< F*(/?x 2 ) = -£- — V V g^. (22) 



(1 + 2k) 2 (-k+ (3^Jk (-1 + k)) K 2 



Consequently, xF*{x) < for x 2 < x < x\. Moreover, G(x) = x 2 /2 and then 
G(x 2 ) < G(xi). Besides, for > 1 we have 



K 2 



(-1 + k) ((-l + 36 2 )K + (e + 6 3 )^/ K (-l + /t))jfi 

< F* (tox) = * J . 

(1 + 2 k) 3 (k + 6y/ K (-1 + k) J 

In other words, F*' (x) > for x\ < x < b\. From Eq. (8) we can see that F*(x) / 
for <| x |-C 1. In addition, since </>'(y) = Cifj, no matter how | y | increases, 
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4>'(y) remains constant. Finally, we can write ^>(x) = F*(x)f(x)/g(x) = (A + B)/C, 
where 



A = 2 + n (l6n 2 + K 2 3 x 3 ) -k (1 + 2k)(xK 2 (3K 2 x + 8k + 4)-12) (23) 
B = 2 In (1 + K 2 x + 2 k) (1 + 2 nf (1 + K 2 x + 2 /c) (24) 



C = s 7r~ o • (25) 

2(1 + 2k) 6 ^ 2 2 (1 + K 2 2; + 2k) 3 

Performing the derivative of ^(x), evaluating it at x = x\a for a > 1, and plotting 
^'(xia), where 



(1 + 2k) 5 (k + o (-1 + k)Y K 2 



^7X^ " - ( 26 ) 



we get the positive function displayed in Fig. 2. Thus, F*(x)f(x)/g(x) is nonde- 
creasing for x\ < x < bi. In this way, we checked each of the three conditions 
of Theorem 1. The unique limit cycle can be seen in Fig. 3 that illustrates the 
numerical simulation corresponding to the results obtained until now. 



6 Case Study 



With the aim of illustrate the results given in the previous sections, we perform 
numerical simulation using the values of the parameters given by Chidambaram 
[5]. All our calculations are performed for a flow characterized by the value of 
the Damkbhlcr number (Da = K\V/F) equal to 300, as reported by Sree and 
Chidambaram [4]. 

Table 2 

The values of the parameters of Cholette's model [4]. 



Symbol 


Value 


Units 


F 


3.333 x 10~ 5 


[m 3 /s] 


V 


io- 3 


[m 3 ] 


K x 


10 


[1/*] 


K 2 


10 


[m 3 /Kmol] 


n 


0.75 


[dimensionless] 


m 


0.75 


[dimensionless] 
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kappa 

sifbn 



Fig. 2. The function flty'^axi), for a > 1, and k > 1. We see that this is a strictly 
positive function. 



Fig. 3. The phase portrait of the Pi-controlled Cholette model subjected to the 
uniqueness conditions. The employed values of the parameters are those given in 
Table 1. 

In Fig. 4 a plot of the evolution in time of the variable x is given for two values of the 
integral parameter Ki . It indicates that the oscillation frequency can be manipulated 
through this control parameter. 




L5 
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Fig. 4. Time evolution of the variable x. The left hand side plot corresponds to 
Ki = 0.5, while the right hand side plot corresponds to Ki = 10. This is a graphical 
representation of the fact that oscillation frequency is a function of the control gain 
parameter Ki. 



7 A Local Center of the Generalized Lienard System ? 



We begin this section by recalling that a limit cycle is an isolated closed 
orbit, while a critical point is a center if all orbits in its neighborhood are 
closed. To the best of our knowledge the literature on period annuli for Lienard 
systems is well developed only for polynomial cases and moreover it focuses 
on Hamiltonian type systems [11], [15], [16]. 
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Notice that for K c (Ref*) = K^(Ref*) + e and e in the following interval 
e* < e < the existence of limit cycles is proved but without guaranteeing 
uniqueness. It is precisely in this interval where our numerical simulations 
point to the existence of a local center in a neighborhood of the origin. Fig. 5 
shows the phase portrait of the Pi-controlled generalized Lienard system with 
a value of K c in the same interval and close to Kf. 
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Fig. 5. The phase portrait of the Pi-controlled generalized Lienard system with a 
value of K^(Ref*) + e in the interval e* < e < and close to Kf, The upper 
plot shows the state space with the limit cycle and the possible annulus in the 
neighborhood of the origin, while the lower plot displays the configuration of the 
vector field close to the origin. 
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8 Conclusions 



The main result of this paper is that the Cholette CSTR model under PI 
control can be mapped into a generalized Lienard dynamical system of non- 
polynomial type. Thus, we establish a new important application of this class 
of nonlinear oscillators that allows us to make a detailed study of the oscilla- 
tory dynamical behavior of these interesting bioreactors. 

Sufficient conditions for the existence and uniqueness of limit cycles of this 
generalized Lienard system are stated in this paper together with numerical 
simulations that indicate the possibility of the existence of a local center (pe- 
riod annulus) when the gain proportional parameter K c of the control law is 
close to the value corresponding to the existence condition of limit cy- 
cles. We also notice that the oscillation frequency is a function of the integral 
control gain parameter Ki, a result that could have practical applications. 
We mention that similar results have been obtained by Albarakati, Lloyd, & 
Pearson [17] for the polynomial case. 

Our work also shows that the Lienard representation of dynamical systems 
and its associated results could have a remarkable potential as an effective 
tool in the control theory for the closed-loop dynamical analysis in the plane. 
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